import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
import math
import random
from scipy.integrate import quad, dblquad
import scipy.stats as st
l_bound = 0
r_bound = 1
x,y = sp.symbols('x y')
#Совместная плотность f(x, y)
f = 2*x**2 + 2*y/3
print("Совместная плотность f(x, y)")
display(f)
#Одномерные законы распределения (маргинальные плотности распределения)
fx = sp.integrate(f,(y, l_bound, r_bound))
fy = sp.integrate(f,(x, l_bound, r_bound))
print("Одномерный закон распределения f_x:")
display(fx)
print("Одномерный закон распределения f_x:")
display(fy)
mult = fx * fy
if mult == f:
print("f(x,y) равное f(x)*f(y) => Х и Y независимыми")
else:
print("f(x,y) не равное f(x)*f(y) => Х и Y зависимыми")
#Условные плотности для непрерывных X и Y
f_xy = f / fy
f_yx = f / fx
print("Условная плотность f(x/y):")
display(f_xy)
print("Условная плотность f(y/x):")
display(f_yx)
def fxy(x, y):
return 2 * (x**2 + y / 3)
def fx(x):
return 2 * x**2 + 1 / 3
def fy(y):
return 2 * y / 3 + 2 / 3
f_max = fxy(1,1)
def neuman_method(N, a, b, f_max):
x_, y_ = [], []
for i in range(N):
while True:
x, y = a + random.random() * (b - a), a + random.random() * (b - a)
z = random.random() * f_max
if fxy(x, y)> z:
x_.append(x)
y_.append(y)
break
return x_, y_
x, y = neuman_method(1000, l_bound, r_bound, f_max)
x_ = np.linspace(l_bound, r_bound, 1000)
fx_x = [fx(x_) for x_ in x_]
plt.hist(x, density=True)
plt.plot(x_, fx_x)
plt.show()
y_ = np.linspace(l_bound, r_bound, 1000)
fy_y = [fy(y_) for y_ in y_]
plt.plot(y_, fy_y)
plt.hist(y, density=True)
plt.show()
Метод Неймана является методом, позволяющим получить значения СВ в соответствии с заданным законом распределения. Этот метод является достаточно универсальным он применим для моделирования всех СВ, значения которых не выходят за пределы ограниченного интервала (a,b).
def xfx(x):
return x * fx(x)
def yfy(y):
return y * fy(y)
def x2fx(x):
return x ** 2 * fx(x)
def y2fy(y):
return y ** 2 * fy(y)
def xyfxy(x,y):
return x * y * fxy(x, y)
Mx = quad(xfx, l_bound, r_bound)[0]
print('Теор. M[x]', Mx)
My = quad(yfy, l_bound, r_bound)[0]
print('Теор. M[y]', My)
Dx = quad(x2fx, l_bound, r_bound)[0] - Mx ** 2
Dy = quad(y2fy, l_bound, r_bound)[0] - My ** 2
print('Теор. D[x]', Dx)
print('Теор. D[y]', Dy)
Mxy = dblquad(xyfxy, l_bound, r_bound, lambda x: l_bound, lambda x: r_bound)[0]
Rxy = (Mxy - Mx * My) / math.sqrt(Dx * Dy)
print('Теор. корреляция', Rxy)
Альтернативная проверка независимости составляющих - так как коэффициент корреляции не равен нулю, то составляющие X и Y являются зависимыми.